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Abstract 



The LSDA+C/ approach to density functional theory is carefully reanalyzed. 
Its possible link to single-particle Green's function theory is occasionally dis- 
cussed. A simple and elegant derivation of the important sum rules for the 
on-site interaction matrix elements linking them to the values of U and J is 
presented. All necessary expressions for an implementation of LSDA+[/ into a 
non-orthogonal basis solver for the Kohn-Sham equations are given, and imple- 
mentation into the FPLO solver ^3] is made. Results of application to several 
planar cuprate structures are reported in detail and conclusions on the interpre- 
tation of the physics of the electronic structure of the cuprates are drawn. 
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1 Introduction 



Density functional theory (DFT) based on the variational principle by Hohenberg 
and Kohn 1 has nowadays a rigorous mathematical basis, mainly due to work 
by Lieb 0. As a theory for (chosen) ground state properties of a many-particle 
system, it holds true for any Coulomb quantum system with arbitrarily strong 
correlations in the ground state. All hampering representability problems of the 
early time are gone (see for instance 0|). However, the central quantity, the 
universal density functional is not known, only its existence can be proved, and 
we have no fully systematic access by approximations. Hence, so far (and very 
likely also in future) we have to model it and to probe the models by comparison 
to phenomenology. This situation is not principally different from other many- 
particle approaches where either models of sufficiently simple Hamiltonians are 
used (in quantum field theory) or the wave function is modeled (for instance in 
Hartree-Fock or Gutzwiller approaches). 

The situation is even less satisfactory in solid state theory, if the focus is 
on the excitation spectra instead on the ground state, because in most cases 
the spectrum of the many-body Hamiltonian has no separate physical relevance 
at all except for its formal use in theoretical expressions for the partition func- 
tion. Instead, what is measured are the spectra of various quasistationary excita- 
tions, defined from few-particle Green's functions, the self-energy parts of which, 
besides being energy dependent non-linear integral operators, are also density 
functionals. 

Although this is not a principal restriction of DFT, the models in use so 
far (local (spin) density approximation, L(S)DA, in the following the acronym 
LSDA is used for both LDA and LSDA, generalized gradient approximation, 
GGA, LSDA plus self-interaction correction, SIC, LSDA plus onsite Coulomb 
repulsion, LSDA+C/, . . . ) are subject to the adiabatic approximation for the 
electron-lattice interaction. 

DFT in the Kohn-Sham (KS) approach to solids yields a KS band structure, 
which as such does also not have a direct physical meaning. Instead, the quasi- 
particle band structure of Bloch electrons is obtained from the self-energy of 
the electron Green's function. It has become common use to speak of weak 
correlations, if in the vicinity of the Fermi level the LSDA KS potential and the 
electron self-energy are not very different. This does by no means imply that the 
correlation energy itself, defined as the difference between the true total energy 
and the Hartree-Fock energy, is small or much smaller than in strongly correlated 
systems (where the LSDA KS potential differs strongly from the electron self- 
energy, the latter often jumps as a function of energy at the Fermi level). It is 
a general experience that in the latter cases also the ground state properties, or 
certain ground state properties, are much worse reproduced by the LSDA than 
in weakly correlated cases (in the above definition). 

This paper deals with strongly correlated systems treated by means of the 
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LSDA+C/ approach. Maybe the first precursor of an LSDA+[/ calculation (at 
that time not fully self-consistent) was the treatment of 4/-electrons by Herbst 
et al. pE|. A first fully self-consistent calculation of values of the Hubbard U was 
performed by Dederichs et al. for the 4/-states of Ce by applying the Koringa- 
Kohn-Rostoker solver of the KS equations to a constraint impurity problem. An 
early similar approach to U for NiO by Norman and Freeman [5] used the APW 
solver in a super-cell formulation. For the further development of this subject 
see [7j. 

Pickett and Wang [H] and Hybertsen and Louie [Hj based the so-called GW- 
approximation for the electron self-energy of semi-conductors and insulators on 
LDA results for the density and KS bands as a starting approximation for the 
GW-approach. In these cases the main difference between the LDA KS potential 
and the self-energy is a jump of the latter, constant in r-space at the Fermi level 
(scissors operation). This is still considered a weakly correlated case, but these 
were the first estimates of the self-energy of an inhomogeneous situation as a 
density functional. 

A systematic incorporation of the Hubbard [/-potential into the DFT model 
functionals started with two papers by Anisimov at al. |10| 111) . While it was 
proposed in ^U] to model the total spin-dependence by the [/-functional and to 
treat the spin independent functional by the LDA, in |llj the orbital polarization 
part (m-dependent occupation of local orbitals) was treated by the [/-functional 
and the isotropic (in r-space) part of the spin-density was treated in LSDA. This 
has the advantage that spin polarization effects can be treated more generally, 
not only in the strongly correlated orbitals. However, this version, later on 
called 'around the mean field', AMF, by Czyzyk and Sawatzky ^2]j gives nearly 
nothing for a half-filled fully spin polarized shell as in Mn 2+ or in Gd. Therefore, 
aiming mainly at reproducing the photoemission spectra (which essentially means 
modeling the electron self-energy rather that the KS potential), an alternative 
[/-functional was introduced in and called the 'atomic limit' version, AL. 
This version which roughly shifts unoccupied orbital energies upward by U/2 
and occupied orbital energies downward by U/2 independent of the shell filling 
(even for filled and empty shells), has been widely used since. For a survey see 

It has to be confessed that all LSDA+[/ models up to now depend on the 
basis used for the KS solver. Most results are obtained so far with an LMTO 
implementation. For a recent APW implementation see |14j . 

In this paper, in Section 2, an FPLO implementation is described. FPLO 
|15j is a high precision high efficiency KS solver which uses a minimum basis 
(and hence is fast) containing only local basis functions which are optimized 
in both a numerical and chemical sense. (It competes in accuracy with well 
converged full-potential APW.) After a short outline how U is integrated into 
DFT, the correlated orbitals used in the FPLO implementation are introduced. 
Since literature statements |12j say that the most important sum rules for the 



4 



interaction matrix elements (screened Slater integrals) of those orbitals are cum- 
bersome to verify, a very simple and elegant derivation is given here. After the 
necessary analysis of the orbital occupation matrix for the non-orthogonal basis 
of FPLO and the introduction of the AMF and AL functionals in the FPLO 
implementation, explicit expressions for the [/-potential and for the total energy 
are given as they are coded in FPLO LSDA+C7. In Section 3, new applications to 
cuprate structures, the 'infinite layer' compound CaCu02, the undoped single- 
layer compound Sr2Cu02Cl2 and the bilayer high-temperature superconductor 
Bi2Sr2CaCu208 are considered and the results are compared to both magnetic 
ground state properties and photoemission spectra. A short summary is given 
in Section 4. 



2 The FPLO implementation of the LSDA+[7 
approach 

The underlying frame of the LSDA+J7 approach is the Hohenberg-Kohn varia- 
tional principle, 

E[v, N] = mm j#[n] + Y, J d 3 rv ss ,(r)n s , s (r) £ J d\n ss {r) = TV j (1) 

for the ground state energy E and spin density ft = (n ss /) of iV electrons in an 
external spin dependent potential v, 



J d 3 rv ss ,(r)n s r s (r) = J d 3 r(vn — B • m), 



(2) 



which holds true in any case of arbitrarily strong correlation. It is based on 
many-particle quantum theory by rigorous mathematics. Of course, the density 
functional H[n] is unknown. 

The generalized Kohn-Sham modeling of this functional is by parameteriz- 
ing the variational spin density by new variational parameters: the Kohn-Sham 
orbitals <fii(r,s) and orbital occupation numbers nj, 

n ss '(r) =^2<pi(rs)ni<p*(rs'), (4>i\4>j) = < n { < 1, ^n i = iV, (3) 

i i 

and by splitting the density functional into an orbital variation expression K[n] 
and a (possibly generalized by gradient terms) local density expression L[n]: 

H[n] = K[n] + L[n] 

K[n] = min I k[<f>i, fanicp* = n \ 

{<f>i,ni} ^ l / -r J > 



(4) 



L[n] = J d 3 rn(r)l(n ss ,(r),Vn,...). 
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This puts the Hohenberg-Kohn variational principle into the Kohn-Sham form 
E[v,N] = mini k[fa,rii] +L Y^n^* + V] rii(fa\v\fa) 

(fa\^ j ) = S ij , o<m<i, J2 n i = N \- ( 5 ) 

While for the L-functional LSDA or GGA models are in use, k is modeled by 
LSDA or LSDA+SIC or LSDA+17. Variation of <p* yields the generalized Kohn- 
Sham equation, 

fa = <pi£i, vf s , = , (6) 

on s i s 

and variation of the ni yields the common aufbau principle which holds true for 
all model variants within this frame and which says that the ground state density 
is obtained by occupying the ./V orbitals with the lowest £j. 
The variants of the LSDA+^7 model correspond to 

k = t + e R + e u , t + e n = ^n i {fa\i\fa) + -^n i n j {fafa\r-^\fafa), (7) 

i ij 

where e u is expressed through projection onto correlated local orbitals \Rfia) 
centered at site (or in the unit cell) R and with orbital and spin quantum numbers 
/j, and a. The projection is given by a local orbital occupation number h^ a which 
depends on the variational quantities fa, nf. 

e=e(n -^' ni]) ' ^ eU = ^dn-^W- (8) 

The functional derivative on the r.h.s. of the last expression yields the projection 
while the partial derivative defines the orbital and spin dependent [/-potential: 
= de u /dn^a-. It is crucial for fitting the models in use into the general 
Hohenberg-Kohn-Sham frame that the correlated orbitals themselves as well as 
the actual value of U are understood fixed and not variational although they 
may be context dependent. They kind of define a location, relevant in a given 
context, in the variational functional space and a functional contribution from 
that location. 

2.1 Correlated orbitals 

In the so called rotationally invariant LSDA+J7 approach the correlated local 
orbitals are assumed to be angular momentum eigenstates centered at R, say, 



1 5k r 
— 777 + * + « 

Hi 0<p* 
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with predefined orbital and spin momentum quantization axes (which both need 
not be the same) 

|.Rjmj0-j), rrn = -k, . . . ,k, (Ti=t,|- (9) 

Only one-site matrix elements, Ri = R2 = R3 = R4, are considered: 

(mim2\w\m3m,4) , w « w(\r — r'|), = 03, 02 = 04 (10) 

The rotational invariance refers to the screened electron-electron interaction, w, 
which is of course an approximation as regards the screening. As a consequence, 
the SO3 transformation properties of the matrix elements are 

( mi m 2 \w\m 3 m^ = ^^(0)^^(0)* 

m' 1 m.2 ,n 3 m 4 

* (m[m 2 \w\m' 3 m'4)U m , 3m3 (d)U m , 4m4 (6), (11) 

where O is any rotation of the r-space and the [/-matrices (not to be confused 
with the Coulomb integral U) yield the relevant SO3 representation: 



/ 



1^(0)17(0) = 1 = u(d)u\d), 

l_, , (12) 
21 + 1 



dO U mim2 {0)U m3Tn4 (O) — S mim4 S m2m ... 



In the last orthogonality relation, dO is Haar's measure of the SO3, J dO = 1. 

These fundamental representation properties allow for a very simple and ele- 
gant derivation of the important sum rules for the matrix elements: Use unitarity 
of U and integrate over dO to obtain 

' S ^{mim2\w\mimi) = 



mi 



£ £ U Lm> \( d ) U Lmt P) {m' im ' 2 \w\m'^ 



m i m^m^m^m'^ 



m2m 2 

X] ^ m ^i 2m ;(0)(^W2l*k3 m 4)^m 4 (0) 



3 Ly m 2 m 2 V^A"''l"''2l" y l"'-3"''4; L/ m^m4V^; ~ (13) 
X] (^^H™! 77 ^) = S m2 m 4 (2l + 1) £A 



The last equation is the definition of the Coulomb integral U. In the same 
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manner, 



mi 



s ^ ( 14 ) 



2/ + 



j ^2 (m^m^wlrn^rn'x) = 6 m2m3 (U + 2U) 



is obtained which additionally defines the exchange integral J. The first result 
(|13|) is intuitively obvious: after summation over mi and integration over r in 
the matrix element, no angular dependence with respect to r' is left except the 
orthogonality (m^m^) = 5 m2 , m4 . The second result (fTl|) is less obvious but 
nevertheless true. 

Expansion of the interaction function into spherical harmonics, 



r i — r 2\) = w ( {r\ + r\ — 2rir2COs9) 



1/2 



A I 



= J2w l (r 1 ,r 2 )Pi{cos9) = ^wi{n,r 2 )— — ^ ^ m (ri)Y^(r 2 ) 

Z=0 «=0 m=-Z 

leads to Slater's analysis 

2li 

(mim2\w\rri3rn4 : ) = .FJ a/ (77117712^37714), 

1=0 

oo 

F t = J J dr x dr 2 {r\Ri(ri)) 2 (r 2 Ri(r 2 )) 2 wi(ri, r 2 ) 
o 

oo 



(15) 



(16) 



y ^ dndr 2 (rxRiin)) 2 ^^)) 2 ^- for / > 0, 



o 



47T 



ai(mim 2 m 3 m A ) = — - - ^ l^ml^msX^imJ^ml^imJ*- 

m=— i 

Here, is the angular momentum of the considered shell, and the second line for 
F[ holds for the unscreened Coulomb interaction which for I > is a reasonable 
approximation since intraatomic screening is effective only for the s-component 
of the interaction. 

Now, from £ m lWr)Y^(r) = P Z (1)(2Z + l)/4vr and 

S ^2ai(mim 2 mim 2 ) = 
4tt 



21 + 1 



^ ] O^hmi | YlO l^krni, 



_ mi 



(^ i m 2 l^ol^ i m 2 )* 



21-4-1 

47T ^ITT 510 (V^ol^mJ* = (2Z * + 1)610 
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it follows immediately that 

U = F . (17) 

Furthermore, 

ai{m l m 2 m 2 mi) = ^ (Yi imi \Yi m \Y hrn2 )(Yi im2 \Y lm \Yi imi ) = 

mim 2 m 1 m 2 



4tt 



mi 



21 + 

Y hm2 {rx)n im2 (r 2 )) (Y, Y lrn{ri)Y? m {r 2 ) 

m 2 m 

(2 |^ )2 1)2 J J dn 1 dn 2 [P li (cos9 12 )] 2 P l (co S 9 12 ) 
(2/j + l) 2 f^ {jD , as , a 



47T 



<*. + !)' I J o „ 



dn[P h (cos9)] Pi{cosO) 



and hence 



(mim2|w|m 2 mi) = (2/j + l) 2 F ; 



mim.2 1=0 



o o oy (is) 

(2/ 4 + l)(C/ + 2/ 4 J). 



Eqs. (|T7|) and (fTH|) relate the Coulomb and exchange integrals U and J to Slater's 
(screened) integrals Fi. 

Recall, however, that the whole analysis presupposes the isotropy of screening 
which could be questioned at least in cases of strong directional covalency. 

2.2 The orbital occupation matrix 

The variants of the LSDA+C/ model are all depending on the basis set of the 
solver of the Kohn-Sham equations. There are a few subtleties in this game which 
never have been discussed in the literature. Here, a non-orthogonal local basis 
implementation ^S] of the solver will be used, since a local orbital representation 
is mandatory for considering strong correlations. Non-orthogonality of the basis 
is rather the rule than the exception for high precision solvers. For an LMTO 
solver see ^H], for an LAPW solver see [14] . 

Consider Kohn-Sham orbitals \k) = \<pk) and orbital occupation numbers 
rik as previously; they need not be eigenstates of spin. Let {|Z)} be a possibly 
non-orthogonal basis for Kohn-Sham orbitals: \k) = J2i \l) c ik, Su> = (Sys- 
tematically, brackets are used for the Kohn-Sham orbitals and parentheses for 
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the local basis orbitals.) For an orthogonal projection onto those basis orbitals 
the contragredient basis \l} = Yli' 10 ' )l'h 010 = is needed. With its 
help, the occupation matrix h = h[4>k,nk] of correlated orbitals \ma) at site R 
in an orthogonal form is introduced as 

flmm'a = ^ y^" 1 ) (Ftma),l (J> I k)n k (k\ l') (S" 1 )z',(Km/cr) = 

k IV ^ 

= ^2 C (Rrna),k n k <?(Bm' cr),k' 
k 

As usually it is assumed that the spin dependence can be made site diagonal by 
choosing a suitable spin quantization axis. The orbital occupation matrix may 
be diagonalized with respect to m, m! at each lattice site R and for each spin 
value a independently: 



» , .%m'a — ^ti ^pLa U ) ■ (20) 



Averages over a correlated shell of angular momentum / 



n °=2l 



^J^fW, n = -(n T + n;), (21) 



are used later on. 

The projector in (jSJ is now 

= E ioc^j^w^ 1 Wytfi*) (22) 

Naturally, in applications the correlated orbitals are assumed to form a subset of 
the basis orbitals, although this is not mandatory. In the FPLO scheme, the basis 
is adjusted in the course of iterations for solving the non- linear Kohn-Sham equa- 
tions. This does not mean that the basis itself is treated as variational. Rather 
the relevant sector of the variational space is tracked along the way of search for 
the Kohn-Sham minimum. Likewise, the relevant location of correlation, that is 
the correlated orbitals as part of the basis, is tracked along. 

2.3 The Orbital Polarization LSDA+£7 Functional 

This functional was introduced under the name 'around the mean field' (AMF) 
in Ref. ^2]. It is zero if the orbitals of an atomic shell are equally occupied, 
hence it depends on orbital polarization. It is given by 



l( n ss'( r ), ■■■) = ^LSDA, 
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e c/,AMF = J ^ a ^_ a \w\^ a ^_ a ){n^ a - n^in^-a - n_ CT ) + 

+ - (Mo-^I^I^Mo-)] (»/i<r - n a )(n^ a - n a ) \ 

1 ■ J]{?7(JV - n CT ) - j(N - n CT ) W, 



(23) 



2 

Ujt 



N a = Y j n^ = {2l + l)fi (T . 



Here, Af is the number of electrons occupying a whole correlated Z-shell, N a is 
that for one spin sort. There is no danger of confusing it with the total electron 
number in Eqs. (^SJ), the latter does not appear any more in the sequel. In the 
second equality use of the sum rules for the matrix elements was made. 

The corresponding fJ-potential is, again most easily with use of the sum rules, 



8e u < 



AMF 



+ - (fi a fi' a \w\n' a fi a )] {h^a - h a ) |> (24) 



One weak point of this version is that it yields no contribution at all in case 
of orbital independent occupation numbers L/„ = n a . This is for instance the 
case of a half- filled completely spin polarized shell (e.g. 4/-shell of Gd). In the 
Gd case this is at least not too bad, as the LSDA gives nearly the right spin 
polarization energy of Gd, although there is a problem with the right magnetic 
ground state (obtained antiferromagnetic in LSDA). 



2.4 The 'Atomic Limit' LSDA+[/ Functional 

With the Gd case in mind and aiming at a better description of the photoelec- 
tron spectra, Czyzyk and Sawatzky introduced another functional in which 
they labeled 'atomic limit' (AL). At least regarding its relation to photoemis- 
sion it should rather be considered a model for the quasiparticle self-energy S 
instead of being related to e . Nevertheless it was given in an e^-form as (again 
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K n ss'(r), ■■■) = ^lsda) 



e U,AL 



\ E j^^-^H/V/O n fMT n lx >- rT + 

^AMF + 1^ ([/ _ J)(1 _ W _ 



(25) 
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The corresponding [/-potential is 

Q e U,AL Q e U,AMF 



(U-J)(h a -±). (26) 



One characteristic feature of this [/-potential is that in case of an isolated shell 
it moves the occupied states downward by (U — J)/2 and the unoccupied states 
upward by (U — J)/2 independent of the shell occupation. By way of contrast, 
the center of the AMF spin subshell potential split moves up with increasing 
subshell occupation (so that the shift of the occupied levels is zero in the case of 
a filled spin subshell and likewise the shift of the unoccupied levels of an empty 
spin subshell; this way yielding no shift at all in Gd.) On the other hand, the 
AL 4/-level splitting of Gd is approximately doubled compared to LSDA which 
is a rather good result in the sense of a self-energy correction. 



2.5 The Kohn-Sham Hamiltonian matrix element and 
the total energy 

Here the formulas are presented, which actually are implemented. For the sake 
of simplicity, the correlated orbitals are identified with selected local basis or- 
bitals as discussed at the beginning of subsection 2.2. The diagonalization of the 
occupation matrix, although greatly simplifying the analytical derivations, does 
not have advantages when coded, since anyway the diagonalizing transformation 
is site- and spin-dependent. Therefore, the full occupation matrix n mm i a is kept. 
(In the following, the site index of all quantities is dropped. If multiple sites 
with correlated states are needed, the formulas apply to all sites separately.) 
The matrix is obtained from the KS states. 

^mrn'a ^ ] c m<r n nka c rn' a (^) 
nk 

Integration is over the irreducible part of the Brillouin zone and a symmetrization 
projector is applied afterwards to get the result for the full zone. 
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The AMF [/-potential matrix as given in the text after Eq. (|SJ) becomes 
v mW = X]X^ f W ,T ' ~ ^"W) [{m' n'\w\mn) - 5 aa r(m'iJ\w\iJ,m)] (28) 

a' fifj,' 

with the property Tr v^^F CT = for each a separately. There is no such property 
of the AL potential since here the up and down shifts are independent of the shell 
occupation. The interaction matrix elements are taken to be spin independent. 
They are calculated from the Slater parameters according to Eq. ()16|) . (Recall 
that the Fi are external parameters, not variational.) 

The projection part of Eq. (J5J) gives only Kronecker deltas. Thus the matrix 
elements of the KS equation are modified by v^^ a for every block of correlated 
orbitals. From the eigenvalues the band structure energy E B is calculated as 
usual. The kinetic energy is obtained from it by subtracting double counting 
corrections. Besides the LSDA corrections the band energy contains a term 

- 2^ V m'ma n mm'a- (29) 
mm' a 

Since the [/-potential matrix is traceless, a constant diagonal term may be added 
to the occupation matrix to obtain 

= } J v m , ma {n mm > a - n ff O mm 'j — IE . {60) 

mm' a 

Finally, the LSDA+[/, kinetic and total energy are 



•C/,AMF _ 


1 ^ J.AMF 




2 


T = 


J.LSDA _ 2£C/AMF 


E = 


^.LSDA _ £>l/,AMF 



The U -potential matrix becomes 

u mW = [(m' fi'\w\mfi) - 5 a(7 i(mfi'\w\fim)] - v^b mm < (31) 

a' fifi' 

with 

v f = U{N-\)-J{N a - 1 -). (32) 

Again, every block belonging to a correlated orbital in the Hamilton matrix is 
modified by adding u^mV- 

Using a nutritive zero it can be shown that 

AT AL ee £ v% ma n mm , a = 2E U ' AL - ^-N (33) 

mm' a 
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Thus, the final result is 



T 



^LSDA _ 2E U ' AL + 




U — J 
2 

U-J, 



N 



E 



^.LSDA _ ^C/,AL 



+ 



N. 
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3 Applications to Cuprates 



The structure of the 2D cuprates considered here may be described by the formula 
B(Cu02)nCa n _i, where B denotes the block layer which separates stacks of n 
Cu02-planes with Ca planes sandwiched in between. 

The LSDA+C/ approach was applied to three compounds (CaCu02 (n = oo), 
S^CuO^Ck} (n = 1) and Bi2Sr2CaCu208 (n = 2)) with the focus on the orbital 
analysis of relevant bands. While the bandstructure gap between occupied and 
unoccupied bands is a ground state property (jump of the chemical potential as 
function of the particle number at zero temperature), the whole bandstructure 
refers to the excitation spectrum, and the KS bandstructure need not compare 
to photoemission, say. Nevertheless, although there is no deeper reason that the 
best LSDA+C/ potential for the KS equations should be close to the electron self- 
energy, there is some hope that like for weakly correlated systems the LSDA+J7 
bandstructure could again provide also an approximation to the quasi-particle 
spectrum. One should, however, be aware that in principle the best KS [/-value 
need not be the same as the best [/-value in Hubbard-type model Hamiltoni- 
ans. The latter value should for instance be used in dynamical mean-field theory 
which is an approach to the electron Green's function (self-energy). A satisfac- 
tory link of the [/-value to the KS variational quantities, the KS orbitals and KS 
orbital occupation numbers, is still missing. For LSDA+[/ as a Hohenberg-Kohn 
model the most relevant results are ground state properties as structural param- 
eters (lattice constants, Wyckoff parameters), magnetic structure and magnetic 
polarization energy and band gap. The lattice parameters for the cuprates are 
obtained in the usual ±2 p.c. agreement with experiment. They will not be 
considered in the following. 

The FPLO version used is 3.00-5 Here, the program settings are summa- 
rized, which are unique to all calculations. The LSDA version is that of Perdew 
and Zunger JB]- The Cu 3d orbitals are taken to be the correlated orbitals. For 
the sake of comparison, the same Slater parameters were used for all calculations: 
U = 8.16 eV (0.3 Hartree), J = 1 eV (F 2 = 9 eV, F 4 = 5 eV). The orbitals were 
optimized in the non-magnetic structures and the resulting compression radii 
were used also for the antiferromagnetic (AFM) LSDA-CZ calculations. Unless 
explicitly else stated the "around mean field" (AMF) functional is used. 
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3.1 CaCuQ 2 



The 'infinite layer' cuprate CaCu02 (no block layer B) does not exist in na- 
ture, but there is an isostructural compound Cao.ssSro.isCuOa- It is considered 
here first because the bands of this infinite cuprate stack are not perturbed by 
hybridization with block layer states and hence in this sense are pure. The 
experimental lattice parameters of Cao.ssSro.isCuC^ JH] are taken for the ficti- 
tious CaCu02- The antiferromagnetic unit cell is shown in Fig. ^ The space 
group is I4/mmm (139). The distance along Cu-O-Cu is d a = 3.86 A and the 
distance in ^-direction between adjacent CuC>2 layers is d c = 3.20 A. (For all com- 
pounds considered here d a is the lattice constant of the non-magnetic cuprate 
plane, while d c is the distance in z-direction between adjacent cuprate multilayers 
[(Cu02) n Ca n _i]. For CaCu02 it is the Cu02 plane distance. It also gives the pe- 
riodicity of the multilayers in z-direction, ignoring a centering shift perpendicular 
to the z-direction.) 

The AFM lattice constants are ao = 6o = V2d a , cq = 2d c . The atom positions 
are Ca (0, \, \), Cu (0,0,0), Cu (0,0, |) and O (f,|,0). Table □ gives the FPLO 
basis set. The number of Fourier components was 1024 per atom and the k-mesh 
subdivision was (6,6,6). Nonrelativistic calculations were performed. 

The symmetry points of the band structures presented below refer to the 
Brillouin zone of the antiferromagnetic cuprate plane (a square of edge length 
2ir/{y/2d a )). For better comparison to the literature we relate it to the non- 
magnetic Brillouin zone which is a rotated by 45° square of edge length 2ir/d a . In 
units of ir /da we have T, (Z) = (0, 0, Q, X, (R) = (±, ± Q and M, (A) = (1, 0, (). 
The first labels refer to £ = and the labels in parentheses refer to £ = 1 = 7r/d c . 
(This is the Z point of a simple tetragonal cell of lattice constant c = d c .) 

Calculations within LSDA and LSDA+C7 were performed. As usual for the 
cuprates, the LSDA gives a metallic ground state. The KS bandstructure is 
shown in Fig. OJ The spaghetti below the Fermi level (here and in all following 
figures put equal to zero) consists of hybridized Cu-3d and 0-2p states, with 
the bonding combinations at the bottom and the antibonding bands formed of 
the orbitals of Fig. [21 crossing the Fermi level. Non-bonding combinations are 
in between. The unoccupied bands above the Fermi level start with Cu-4s and 
Cu-4p character and then enter a bunch of Ca-3d bands above 5 eV. 

The same bands weighted (by linewidth) with the square of the coefficient of 
selected basis orbitals in the KS state are shown in Fig. |1J From the upper two 
panels one can read off a (Cu-O)^ covalency split of more than 5 eV while the 
third panel shows a (O-O) covalency split of about 3 eV. One further observes 
that the CT orbitals and the O z orbitals hybridize also with the Cu 4s and 4p 
orbitals while the in-plane O^ orbitals hybridize additionally with the Ca 3d 
orbitals. Figs. El and 0] are presented here for comparison with the LSDA+J7 
results shown below. 

Experimentally, (Cao.85Sro.i5)Cu02 is an AFM insulator with a bandgap of 
more than 1 eV and a Neel temperature T\r ~ 540 K. The [/-functional cures 
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this deficiency and one finds an AFM solution with a spin polarization energy 
AE = 27.9mHa per formula unit below the Pauli-paramagnetic (PM) state. The 
site projected copper 3d moment is 0.71/xb and the total copper spin moment is 
0.69 fjis (reduced by negative 3s3p moments). 

The already discussed two relevant molecular orbitals (MO) of Fig. El are the 
candidates for the highest occupied molecular orbital (HOMO) in the correlated 
electronic structure of cuprates, that is, those MOs which are relevant for the 
valence band edge. The LSDA+C/ KS bands of CaCu0 2 with the AMF and AL 
functionals are shown in Fig. El The main difference between both functionals 
is found in the unoccupied bands. As to be expected for a more than half-filled 
shell, the upper Hubbard band of (Cu-3d x 2_ y 2 and 0-2p a character, see below) 
lies higher in the AMF case compared to AL. In AL this band is the lowest 
unoccupied band. In the occupied part near the Fermi level the differences are 
small, while far below larger differences are found due to the different position 
of the lower Hubbard band and thus to different hybridization. Since this paper 
focuses on the occupied bands near the Fermi level, the differences are not very 
relevant, and all further results are presented for the AMF functional. 

On Fig. H3 the orbital weights to the AMF bands in analogy to Fig. 0] is 
shown. First, it is clearly seen that the upper and lower Hubbard band is formed 
by Cu-3d x 2_ y 2 and 0-2p a orbitals. At the valence band edge (points X and R 
shown, but without noticeable dispersion on the whole line X-R in ^-direction of 
the fc-space) 0-2p a and 0-2^ contribute equally strongly while the Cu-3d x 2_ y 2 
orbital contribution is largely suppressed compared to the LSDA result. This 
suppression of the Cu contribution to valence holes in cuprates is confirmed by 
experiment [21]. The 0-2p n contribution on the other hand is strongly enhanced 
compared to LSDA, a new result which is missed in most model Hamiltonian 
treatments in the literature where the 0-2p n degree of freedom is excluded from 
the Hamiltonian in most cases. Other orbitals do not contribute. Hence, on 
the line X-R the highest occupied band is mainly a hybrid of 0-2p a and O^p^. 
These band states are assumed to form the Zhang-Rice singlet with the nominal 
3<i-hole on the Cu site (upper Hubbard state) |21j . The bandwidth of the highest 
valence band is about 1 eV due to hybridization with a flat band 1 eV below the 
Fermi level. On the line X-M (R-A) the 0-2p a contribution fades away towards 
M (A) due to this hybridization. In LSDA the 0-2^ orbitals do practically 
not contribute to this band which led to the neglect of that orbital in model 
Hamiltonians. This failure of the LSDA is mainly due to the fact that in LSDA 
the 0-2p w bands are deeper in energy compared to the Cu-3(i a ,2_j / 2-0-2j? (7 bands. 

The absolutely highest occupied band in the LSDA+C/ result is a pure 0-2^ 
band (at T) which therefore is in energetic competition with the Zhang-Rice state 
when adding additional holes. The 0-2p n bands show a considerable dispersion 
in z-direction which comes from a weak hybridization with the unoccupied Ca- 
3d x 2_ y 2 orbitals. The Cu 4s orbitals do not contribute at the valence band edge. 
They are mixed into the unoccupied states and are also slightly mixed into 0-2p a 
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states about 1.5 eV below the Fermi level on the line M-A. 

3.2 Sr 2 Cu0 2 Cl 2 

The single layer compound Sr2Cu02Ci2 is probably the most two-dimensional 
of all cuprates. The block layer consists of 2 SrCl layers which separate single 
Cu02 planes. Between these SrCl layers crystals are easily cleaved, whence 
most photoemission data on undoped planar cuprates are recorded from this 
material. Two adjacent cuprate planes are shifted horizontally by a shift vector 
(5,5,0) relative to each other, which produces a body centered non-magnetic 
unit cell. The large distance of cuprate planes from each other prevents valence 
state coupling in z-direction. The AFM cell has base centered orthorhombic 
symmetry and is shown in Fig. 

The standard cell choice for the orthorhombic cell, which has also to be 
used in the FPLO code results in c-base centering. (The stacking direction is 
the 6-direction and the cuprate plane is the c,a plane.) The space group is 
Cmmm (65). The cell parameters are d a = 3.973 A and 2d c = 15.618 A |25j . 
The lattice constants of the AFM cell are ao = cq = V2d a , bo = 2d c . The 
atom positions are Cu (0,0), Cu (0,0,0), O (-|,0,-|), CI (0,-0.317, 5), CI 
(0, 0.183, 0), Sr (0, -0.107, |) and Sr (0, 0.393, 0). TableElgives the basis set. The 
number of Fourier components was 1000 per atom and the k-mesh subdivision 
was (12,12,12). Scalar relativistic calculations are performed. For the sake of 
comparison the stacking direction is further on denoted z and the cuprate plane 
as the x, y-plane. The symmetry points are chosen corresponding to the scheme 
described in the previous subsection. 

The site projected copper 3d moment of the LSDA+C/ result is 0.758/Ub and 
the total copper spin moment is 0.748/i£- In Fig. |H] the orbital weights for the 
relevant bands are shown. The situation is rather similar to that of CaCu02, so 
only the differences are pointed out. 

The "Zhang-Rice" band is slightly higher and touches the Fermi level at X 
and R. The intersecting band is correspondingly lower at the line M-A, as a 
consequence the width of the upper valence band is about 1.5 eV and the fading 
of the 0-2p a character on the lines X-M and R-A is less pronounced, a reduction 
of about 50 p.c. remains. The z-dispersion of the 0-2^ band on the line T-Z has 
gone as there are no Ca-3d states present for hybridization. There is a marked 
2D character of the compound. 

The quasiparticle low energy dispersion measured by ARPES (single hole 
excitation) is shown in Fig. El A detailed discussion is found in 26 and [2*7] . 
The left part of the experimental spectra (from (0,0) to (5, 5)) compares nicely 
with the 0-2p a dominated LSDA+J7 band on the line T — X (second panel of 
Fig.IHJ) and the right part (from (§,§) to (0,7r)) with the same LSDA+C/ band on 
the line X — M. Even the reported fading ARPES intensity when going from X 
towards M (see also [28 ) agrees with the fading 0-2p a projection of that band. 
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Nevertheless, the experimental band width is smaller by a factor of about two to 
three and the situation on the line r — M is less clear although the comparison 
of only the 0-2p a projected bands (second panel of Fig. El cf. the discussion 
of models above) to model results contained in Fig. is not so bad. After all, 
LSDA+C7 accounts for electron correlations still rather grossly. 

3.3 Bi 2 Sr 2 CaCu 2 08 

The bilayer compound Bi 2 Sr 2 CaCu 2 0s has a double layer of cuprate planes sep- 
arated by block layers. The block layer consists of two BiSr0 2 layers between 
which crystals cleave equally easily as in the previous case and also favor the ma- 
terial for photoemission. Adjacent double layers are again shifted horizontally 
by a shift vector (^, ^,0) relative to each other, which produces a body centered 
tetragonal non-magnetic cell. The hypothetic AFM cell again is assumed base 
centered orthorhombic. 

The cell parameters are again given in the c-base centered setting, result- 
ing in the 6-axis being the stacking direction. The space group is Cmmm (65). 
The cell parameters are d a = 3.817 A and 2d c = 30.6 A. The lattice constants 
of the AFM cell are ao = cq = V2d a , &o = 2d c . The atom positions are Ca 
(0, \\), Sr (0, 0.3597, \), Sr (0, 0.1403, ±), Bi (0, -0.4478, \), Bi (0, -0.0522, ±), 
Cu (0, -0.3040, I), Cu (0, -0.1960, \), Oi (-±,0.1960, 0), 2 (0, -0.1250, ±), 2 
(0,-0.3750, \), 3 (0,0.0450, \) and 3 (0,0.4550, \) . The two Sr, Bi, 2 , 3 
and Cu atoms are equivalent in the non-magnetic cell. The Oi atom is that of 
the cuprate plane. The two Cu spins are antiferromagnetically ordered in the 
AFM cell. Table 01 gives the basis set. The number of Fourier components was 
500 per atom and the k-mesh subdivision was (8,8,8). Again, scalar relativis- 
tic calculations were performed and the symmetry points in fc-space are chosen 
corresponding to the scheme described in subsection 3.2. 

For this material both the LSDA and LSDA+[7 yield a metallic ground state 
with Cu-0 bands and Bi bands crossing the Fermi level. While the LSDA re- 
sults in a non- magnetic solution, the LSDA+C7 calculation yields a stable AFM 
state. The site projected copper 3d moment is 0.696/^b and the total copper spin 
moment is 0.684/i#. 

Fig. ITU1 shows both the LSDA and LSDA+C/ bands of Bi 2 Sr 2 CaCu20 8 close 
to the Fermi level. For a better orientation in the band character the Cu-3d c 2_ J/ 2 
projected bands are also shown. Note that the LSDA+J7 results on the right 
panels have twice as many bands as the LSDA results on the left panels due 
to the AFM order of the former ground state. The bands crossing the Fermi 
level and not seen in the lower panels are Bi bands (more precisely BiO bands 
hybridized with orbitals of the block layer oxygen). 

The splitting of the Cu-3dj.2_.y2 projected LSDA bands crossing the Fermi 
level on the line X-M (R-A) is the much discussed bilayer splitting between bond- 
ing and antibonding combinations of the CuO states in both Cu0 2 layers of the 
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bilayer |27] . The coupling of those states is mainly due to a small hybridization 
with Cu-4s states, and this part of the coupling has k^ — ky symmetry. Hence, the 
splitting is maximum (~ 0.25 eV) at point M (A) and nearly zero at point X (R) 
although a very small splitting remains there due to a small direct coupling. It is 
readily seen by a simple symmetry argument that in the AFM state this splitting 
must be zero on the whole line X-M (R-A), if spin-orbit coupling is neglected. 
Accordingly all bands on the right panels of Fig. El] are twofold degenerate on 
these lines and do not show a bilayer splitting there. This only develops away 
from these lines, for instance in the X-r and M~r directions. However, the Cu- 
2>d x 2_ y 2 projected LSDA+[7 bands on the right lower panel show around M and 
A another splitting of approximately the same magnitude which is due to the 
crossing of another oxygen band. This point has never been considered in the 
literature to date. 

Note also that the non-magnetic structure of the left panels has a larger Bril- 
louin zone so that the rising band from T to X continues to rise from X=(^, ^,0) 
to the Brillouin zone corner (1,1,0), whereas the AFM structure of the right pan- 
els has a charge transfer gap (dominated by U and hence much larger than the 
exchange splitting) at X, which now lies on the Brillouin zone boundary, and the 
conduction band has its maximum at X. On the other hand, the band pair close 
to the Fermi level on the line r — Z on the upper right panel is on the line (1,1, C) 
in the non-magnetic state and not shown on the upper left panel. To illustrate 
the charge transfer gap, the LSDA+t/ band structure projected on Gu-3d x i _ y 2 
is shown once more in Fig. for a larger energy window. 

Experimentally, the oxygen in the block layers of Bi2Sr2CaCuOs is volatile 
and its stoichiometry is governed by thermodynamics. Moreover, the geometry 
of the block layer is distorted in a disordered way compared to the ideal structure 
used in the calculation. In the recent, highest resolution photoemission spectra 
|27j . Fermi surface pockets around point M which should be present due to the 
Bi bands in both the LSDA and LSDA+[7 results are not seen. Hence, one could 
assume that they are pushed away from the Fermi level (together with the oxygen 
bands on T-Z) by a distortion potential. Correspondingly, the Fermi level would 
be lowered to ensure the electron count. This would be a big problem for the 
LSDA band structure where not only the antibonding bilayer split band would 
have an electron Fermi surface closed around T but also the bonding one, both 
in contradiction to what is seen in photoemission. The LSDA+J7 band structure 
(lower right panel of Fig. HOD on the contrary would be in rather good agreement 
with photoemission if one wipes out the down folded bands of the AFM order 
which is not found in experiment. (The material cannot be reduced down to 
undoped CuOg - planes.) 

There would be another stark difference to what is discussed in the literature 
with respect to the bilayer splitting: If the band splitting around M and observed 
in photoemission would be a bilayer splitting, one should expect it to be strongly 
reduced when reducing the doping level to an underdoped superconductor which 
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is regarded to develop strong AFM correlations (also seen in neutron scattering). 
It should be reduced to zero where AFM order sets in. The splitting of the 
LSDA+C/ bands of different origin on the contrary is to be expected largely 
independent of antiferromagnetic order and hence on doping. 

3.4 Implications on Magnetic Interactions 

Magnetic couplings may grossly be obtained from total energy differences of 
LSDA+J7 results for ferromagnetic and antiferromagnetic order; in more detail 
they may be obtained from calculated energies of spin spiral states £22 • 

For a more detailed understanding of their physics, a tight-binding model 
of the kind of Emery's model should be extracted from the LSDA+C/ results, 
which then may be down-mapped to a kind of a t — J model. For a hypothetic 
ferromagnetic order (assumed for the sake of simplicity) such a tight-binding 
model was derived in Ref. 23 . 

One main conclusion from the present orbital analysis is that in a large part 
of the Brillouin zone there is a strong hybridization of 0-2p a with 0-2p n orbitals 
in bands hybridized with the Cu-3d,,2 _ y 2 orbital. Furthermore, the 0-2p n bands 
are in energetic competition with the bands forming the Zhang-Rice state, if 
additional holes are doped. Hence, the 0-2p w orbitals must be included in the 
Emery model in order to correctly describe the t-terms which determine the 
magnetic coupling. 

In CaCu02 there is a sizable dispersion of the in-plane 0-2p n bands in z- 
direction mediated by some hybridization with Ca-3d x 2_ y 2 orbitals. This is in 
accordance with the experimental finding that CaCu02 has the highest Neel 
temperature, Tjy = 540 K, of all layered cuprates indicating 3D magnetism |19[ 
1317] . For a survey of the magnetic properties of cuprates see Ref. [3*T] . 

By contrast, in S^CuC^C^ there is no dispersion of the corresponding bands 
in z-direction due to the Sr2Cl2 buffer layers. There is only dipole-dipole coupling 
of the planes, compatible with the experimental findings, T/v = 256 K |32j . 

4 Conclusions 

The LSDA+?7 approach is shown to fit perfectly in the frame of DFT by Hohen- 
berg and Kohn for the electronic ground state, provided the theory can be closed 
by linking the [/-value to the variational quantities, the KS orbitals and orbital 
occupation numbers (which explicit link is yet to be rendered). On the other 
hand, the LSDA+C/ potential is also widely understood as an approximation to 
the electron self-energy (first step towards an LSDA+dynamical mean-field the- 
ory). There is hope that in this way like in the situation of weakly correlated 
systems it provides again a tool to obtain a rather accurate ground state and a 
reasonable approximation to the quasi-particle spectrum (band structure) in a 
single run from only one set of equations. To pursue this goal, in the first part of 
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the present paper an attempt was made to present the structure of the LSDA+J7 
theory as clearly as possible. 

Application to several typical planar cuprate structures and comparison to 
experimental data of the isolating electronic state and the magnetic state as 
well as of quasi-particle spectra probed by photoemission seems to support this 
expectation. On the other hand it revealed a number of new aspects in the physics 
of the electronic structure of cuprates. Notably the sufficient completeness of 
most model Hamiltonians in use must be questioned, at least in connection with 
the dimensionality and details of magnetic couplings, and the so-called bilayer 
splitting of the band structure of E^S^CaC^Os should be reanalyzed. 
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Atom 


core 


valence 


Ca 


Is 2s 2p 


3s(-l) 3p(-l) 4s(1.1240) 4p(1.0927) 3d(-l) 


Cu 


Is 2s 2p 


3s(-l) 3p(-l) 4s(1.3148) 4p(1.2648) 3(2(1.3500) 





Is 


2s(1.2869) 2p(1.2822) 3d(-l) 



Table 1: Basis set for CaCu0 2 . (Compression parameter for valence orbitals in paren- 
theses.) 



Atom 


core 


valence 


Cu 


Is 2s 2p 


3s(-l) 3p(-l) 4s(1.2300) 4p(1.1199) 3(2(1.3148) 





Is 


2s(1.2580) 2p(1.2377) 3(2(1.1162) 


CI 


Is2s2p 


3s(1.1518) 3p(1.1338) 3(2(1.0852) 


Sr 


Is 2s 2p 3s 3p 3d 


4s(-l) 4p(-l) 5s(1.1360) 5p(1.0816) 4(2(1.1503) 



Table 2: Basis set for Sr 2 Cu0 2 Cl2. (Compression parameter for valence orbitals in 
parentheses.) 



At 



om core valence 



Ca Is 2s 2p 3s(-l) 3p(-l) 4s(1.1298) 4p(1.1210) 3(2(1.2286) 

Sr Is 2s 2p 3s 3p 3d 4s(-l) 4p(-l) 5s(1.1239) 5p(1.0717) 4(2(1.1435) 

Bi Is 2s 2p 3s 3p 3d As Ap Ad Af 5s(-l) 5p(-l) 6s(1.3570) 6p(1.2922) 5(2(1.4752) 

Cu Is 2s 2p 3s(-l) 3p(-l) 4s(1.2229) 4p(1.1549) 3(2(1.3060) 

01 Is 2s(1.2464) 2p(1.2177) 3(2(-l) 

2 Is 2s(1.1766) 2p(1.1895) 3(2(-l) 

Q 3 Is 2s(1.0683) 2p(1.0951) 3(2(-l) 

Table 3: Basis set for Bi 2 Sr 2 CaCu 2 8 . (Compression parameter for valence orbitals in 
parentheses.) 
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Figure 1: Unit cell of antiferromagnetic CaCu0 2 ; Cu 2 spin up, Cu 3 spin down. 




Figure 2: Left: Antibonding (k = (n,ir,0)) dpa-orbital commonly assumed as the 
HOMO that forms the Zhang- Rice singlet together with the nominal Cu-d hole |27Hl21j. 
Right: 0-0 antibonding (k = 0) in-plane ]?7r-orbital, lifted up by crystal field and 
weakly hybridized with Ca-d orbitals in adjacent layers: the true HOMO of the 
LSDA+^7 model [221 123 • 



24 



CaCu0 2 I LSDA 
10.0 




txm r zra z 



Figure 3: LSDA band structure of CaCu02- 
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Figure 4: LSDA result for CaCu0 2 . From top: orbital weight of the Cu-Sg^ 2 -?, 2 orbital, 
the 0-2^ orbitals, the 0-2^ and the 0-2p z orbitals. 
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CaCu0 2 | LSDA+U (U=8.16 eV, J=1) M Cu =0.71 
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Figure 5: LSDA+c7 band structure of CaCu0 2 . Top: AMF, bottom: AL. 
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CaCu0 2 | LSDA+U (U=8.16 eV, J=1) M Cu =0.71 
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Figure 6: From top to bottom: orbital weight of the C\i-3d x 2_ y 2 orbital, the 0-2p a 
orbitals, the 0-2p w and the 0-2p z orbitals. 
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Figure 7: Unit cell of antiferromagnetic Sr 2 Cu0 2 Cl2; Cui spin up, Cu 2 spin down. 
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Figure 9: Energy dispersion of quasiparticles for Sr 2 Cu0 2 Cl 2 |26j. The energy zero is 
put at the top of the band, about 0.7 eV below Fermi level. Open symbols: experimen- 
tal data; solid circles: self-consistent Born approximation for a t — t' — t" — J model; 
solid line: tight-binding fit; dashed: t — J model; dotted: spinon model dispersion. 
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Bi 2 Sr 2 CaCu0 8 | LSDA 



Bi 2 Sr 2 CaCu0 8 | LSDA+U (U=8.16 eV, J=1) M Cu =0.696 
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Figure 10: Bi 2 Sr 2 CaCu208: left: LSDA band structure and Cu-34 2 - y 2 orbital weig 
right: LSDA+£/ AMF bandstructure and weights. 
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Figure 11: Bi 2 Sr 2 CaCu 2 8 : LSDA+£/ AMF Cn-M x 2_ y 2 orbital weights. 
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